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Abstract. Many modified gravity theories are under consideration in cosmology as the source 
of the accelerated expansion of the universe and linear perturbation theory, valid on the largest 
scales, has been examined in many of these models. However, smaller non-linear scales offer 
a richer phenomenology with which to constrain modified gravity theories. Here, we consider 
the Hu-Sawicki form of f{R) gravity and apply the post-Friedmann approach to derive the 
leading order equations for non-linear scales, i.e. the equations valid in the Newtonian-like 
regime. We reproduce the standard equations for the scalar field, gravitational slip and the 
modified Poisson equation in a coherent framework. In addition, we derive the equation for 
the leading order correction to the Newtonian regime, the vector potential. We measure this 
vector potential from f{R) N-body simulations at redshift zero and one, for two values of 
the /rq parameter. We find that the vector potential at redshift zero in f{R) gravity can be 
close to 50% larger than in GR on small scales for I/rqI = 1.289 x 10“®, although this is less 
for larger scales, earlier times and smaller values of the fno parameter. Similarly to in GR, 
the small amplitude of this vector potential suggests that the Newtonian approximation is 
highly accurate for f{R) gravity, and also that the non-linear cosmological behaviour of f{R) 
gravity can be completely described by just the scalar potentials and the f{R) field. 
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1 Introduction 

A key problem in modern cosmology is understanding the origin of the late time accelerated 
expansion of the universe. One of the possible solutions to the problem is that General 
Relativity (GR) is not the theory of gravity that operates on larger scales in the universe. 
The accelerated expansion aside, applying GR to the entire universe is an extrapolation over 
many orders of magnitude compared to the scales where it is tested. For these reasons, there 
is increasing interest in considering alternative gravity theories and their phenomenology, see 
[8] for a review. 

On larger scales, where inhomogeneities are small and standard perturbation theory can 
be applied, much work has gone into showing that for many models the extra effects can be 
described by a time- and space-dependent Newton’s constant and “gravitational slip”, i.e. a 
non-zero difference between the scalar potentials, see e.g. [13, 28, 36] or further references and 
discussion in [8, 9]. In addition, analytic forms can often be found for these effects. On smaller 
(non-linear) scales, there is potential for more phenomenology, due to the intrinsically non¬ 
linear nature of gravity, screening mechanisms (see e.g. [16]) such as the chameleon mechanism 
[17, 23] and the possibility of significantly sourcing vector and tensor modes. 

Non-linear scales are typically studied through the use of N-body simulations, either 
Newtonian for the case of GR, or “Newtonian-like” for the case of modified gravity theories 
where the Newtonian gravitational equations are modified. It has been an ongoing problem 
in cosmology to understand how, for a GR universe, the non-linear Newtonian gravitational 
equations arise, how well their solution matches that of a GR cosmology (see e.g. [10]), and 
how we can move beyond the Newtonian simulations to simulations that encapsulate more 
of GR (see e.g. [1, 26]). However, little attention has been given to the same questions in 
modified gravity cosmologies. In addition, the derivation of the “Newtonian-like” equations 
for different modihed gravity theories can be somewhat incoherent, in the sense that several 
different arguments and approximations are invoked in order to simplify the equations and 
remove many of the terms. 
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One of the popular modified gravity models under consideration is f{R) gravity, first 
suggested as a possible cause of the accelerated expansion in [4, 5]. In this model the GR 
Lagrangian, consisting solely of the Ricci scalar R, is modified to have an additional function of 
the Ricci scalar f{R) included. Of particular interest is the Hu-Sawicki form for f{R) gravity 
[12], which has several nice properties: It can be tuned to match a ACDM expansion history 
and can also evade Solar-System tests through the chameleon [17, 23] screening mechanism. 
N-body simulations have been run for the Hu-Sawicki form of f{R) gravity [25, 34] and truly 
non-linear phenomena have been observed, including the manifestation of the chameleon 
screening [34]. As a result of this chameleon mechanism, the f{R) modifications to the GR 
equations disappear in dense environments, which is the reason why f{R) gravity can pass 
Solar-System tests. 

In this paper, we will use the post-Friedmann formalism [21, 22], a post-Newtonian type 
expansion in powers of the speed of light c, designed for a FLRW cosmology. We will apply 
this formalism to f{R) gravity in order to examine the Newtonian regime for this theory. In 
particular, we will derive the leading order equations and the first correction to the Newtonian 
regime, a constraint equation for the vector potential. We will use these equations to calculate 
the vector potential from N-body simulations and use this to comment on the accuracy of the 
“Newtonian-like” simulations for f{R) gravity, similarly to the work carried out in [3, 31] for 
GR. 

This paper is laid out as follows. In section 2 we review the pertinent details of Hu- 
Sawicki f{R) gravity and the post-Friedmann formalism and in section 3 we apply the post- 
Friedmann formalism to f{R) gravity. In section 4 we calculate the vector potential from 
f{R) N-body simulations, compare it to the GR case and comment on the consequences. We 
conclude in section 5. Throughout the paper, a horizontal bar denotes background quantities 
and V^A = A^a, where the indices on the partial derivatives can be raised and lowered freely. 
Dots denote partial derivatives with respect to time and a subscript 0 denotes quantities 
evaluated at redshift 0, i.e. today. Our metric sign convention is (—h 


2 Review of relevant theory 


2.1 f{R) gravity 

The equations of f{R) gravity are derived by generalising the Lagrangian of GR, such that a 
generic function of the Ricci scalar R, f{R), is added to the action, 

c_ f + I c 

S — Id Xy g *5niatter* V^-l) 


Specifying the function f{R) then completely specifies the theory. In general, f{R) gravity 
contains fourth order equations of motion, and is a subset of scalar-tensor theories. The 
modified Einstein equations, derived by varying the new action with respect to the metric 
are given by 


SvrGV = (1 + fR)-^{R + f{R)) + g^uDfR - V^V,/k, (2.2) 

where is the covariant derivative, □ = and /r = df{R)/dR. Jr corresponds to 

the extra scalar degree of freedom that is present in f{R) gravity. The Hu-Sawicki form [12] 
of f{R) gravity is given by 
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The ratio Ci/c 2 is constrained by reqniring that the f{R) theory reprodnces the ACDM 
expansion history, Tr?c\lc 2 = IbvrG/OA- The equations can be re-written in terms of /ijp, the 
valne of the field in the backgronnd today, 

fR = fRo (§) ^ (2.6) 

/(i?) = -167rGpA-^^ (2.7) 

Typically, the free parameter n is taken to have the valne nnity, leaving as the only free 
parameter in the theory. As mentioned in the introdnction. This form of f{R) gravity can 
match a ACDM expansion history and contains the chameleon screening mechanism [17, 23]. 
In dense environments, this mechanism screens the effects of the fifth force in the modified 
Einstein eqnations, either fnlly or partially. This allows the Hn-Sawicki model to fit Solar- 
System constraints [12]. 


2.2 Post-Friedmann formalism 

The post-Friedmann formalism was proposed in [21, 22] and has been nsed to examine the 
vector potential in a GR-I-ACDM cosmology [3, 31] and also the weak-lensing deflection 
angle for non-linear scales [30]. It comprises a post-Newtonian type expansion of the Einstein 
eqnations in powers of the speed of light c, altered compared to a Solar-System type expansion 
in order to apply to a FLRW cosmology. The pertnrbed metric is considered in Poisson gange 
and expanded np to order c~^, 
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The two scalar potentials have each been split into their leading order (Newtonian) {U]sf,VN) 
and post-Friedmann (Up,Vp) components. The vector potential in the Oi part of the metric 
has also been split np into and Bf. As this metric is in the Poisson gange, the three 
vectors B^ and Bf are both divergenceless, B^- = 0 and RT = Q. In addition, the tensor 
pertnrbation hij is transverse and tracefree, h] = hfj = 0. From a post-Friedmann view¬ 
point, the terms of order c~^ and c~^ are considered to be leading order and the terms of 
order c~‘^ and c~^ are considered to be higher order. As is nsnal, the time derivatives will 
also come with a factor of c~^. The formalism was designed for a ACDM cosmology, so the 
energy-momentnm tensor is constrncted from the fonr-velocity of pressnreless dnst and then 
expanded in powers of c. See [22] for details of the energy-momentnm tensor, as well as nsefnl 
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expressions for the Ricci scalar and the other important quantities. In this paper, p, 6 and 
Vi denote the background density, density perturbation and velocity for the pressureless dust 
fluid. 


The leading order equations (c“^ order) derived using this formalism are equivalent 
to the quasi-static, weak field and low velocity limit of the Einstein (or modified Einstein) 
equations. The advantage of the post-Friedmann expansion is that the corrections to this 
limit can be examined order by order. In a GR+ACDM cosmology, the first correction is a 
constraint equation for the vector potential. 
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with It, denoting the vector (divergence-less) part. This was used to calculate the vector po¬ 
tential from ACDM N-body simulations in [3, 31]. This vector potential acts as a quantitative 
check of the Newtonian approximation, and could also influence cosmological observables such 
as weak-lensing [30]. 


3 Post-Priedmann f{R) gravity 


We will now derive the leading order gravitational equations valid for non-linear scales in 
f{R) gravity. These equations come from expanding the modified Einstein equations using 
the post-Friedmann formalism. The Ricci scalar, Ricci tensor, metric and energy-momentum 
tensor will all be expanded as in GR. In addition, we need to decide how to expand the f{R) 
and fji terms. We have some clues as to how to perform this. Since R ~ at leading order, 
we would expect f{R) to be a factor of c~^ smaller than //j at leading order. Examining 
the linear perturbation equations, and those currently used in f{R) N-body simulations, the 
fji field behaves similarly to the scalar potentials in the metric, i.e. the value of the field is 
small everywhere but its spatial derivatives can be large. This suggests that at leading order, 
the fa field should be of similar order to the scalar potentials. Finally, considering the form 
of the f{R) function, we would want the leading order perturbation to occur at the same or 
higher order than the background term, as is the case for the Ricci scalar and Ricci tensor. 
These considerations can all be satisfied by expanding the terms at leading order as follows. 
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In expanding the terms like this, we have assumed that the weak field approximation holds 
for the full field, both the perturbations in this field and its background value. This 
holds for the cosmologically interesting Hu-Sawicki models that we are concerned with here. 
In addition, we have assumed that the quasi-static approximation will hold. This has been 
investigated previously and found to hold on the scales that we are dealing with here [2, 24]. 
We now expand the modified Einstein equations to generate the leading order equations. 
As well as this ansatz, we note that at leading order Rq/R = 0(1) and Rq = 0(c~^). We 
subtract the homogeneous background from the equations and define SR* = R* — R* and 
(5/^ = f'^iR) — f'^{R). Here, SR* and <5/^ are not linearised; they are the full perturbed 
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quantities, with only the background subtracted off. For the covariant derivative, we note 
that □ can be replaced with as the time derivative terms will be of higher order. More 

generally, when two covariant derivatives are acting on a scalar A, VqVqA and ViVjA can 
both be replaced by the partial derivatives, as the additional terms that appear are higher 
order. However, VgViA = A^i — ^A^i, where now the additional term from the covariant 
derivative is at the same order so must be included. 

Then, the trace equation is 

RdfR + 6R{1 + MR) + Sfn) -2{5R + 5f{R)) + SDJ/k = ^<5r. (3.3) 

As we are only interested in the leading order equations here, we then expand each term as 
described earlier to show the leading order power of c for each term 
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We then repeat the same process for the 00 Einstein equation, 

500 
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Expanding out the terms and then taking only the leading order yields 
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This equation can be combined with the trace equation to give 
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which will we use shortly. A further form for the f{R) Poisson equation is given by 

^^^pb = kv^bfR - \MUn. 
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Finally, we need to examine the diagonal ii Einstein equations in the same way (note that 
there is no snmmation here yet), 


RiiSfu + 6Rii (1 + fR{R) + dfR) - ^ (<5i? + 5f{R)) -'^{R + 6R + f{R) + 5f{R)) 
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Again, expanding ont and then keeping only the leading order terms gives 
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We can then snm the three ii eqnations and combine with eqnation 3.8 to give 


4 v' {Un - Vn) = 

The final set of leading order eqnations is thns 
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As expected, these post-Friedmann eqnations match those cnrrently nsed in f{R) N-body 
simnlations. In particnlar, as 5R* and <5/^ are not linearised, the fnll chameleon mechanism 
is incorporated in these eqnations. The effects of the chameleon mechanism have been seen in 
f{R) N-body simnlations that implement these eqnations, see e.g. [35]. Note that in deriving 
eqnations (3.13), we have only expanded ont the Ricci tensor part of the Einstein tensor. The 
pertnrbed Ricci scalar has been kept in the eqnations and then nsed to snbstitnte for the /r 
term. This is why the potential Vn that appears in the Poisson eqnation in [22], here only 
shows np in the last of eqnations (3.13) and not in the Poisson eqnation. For 5/^ = 0, the 
last eqnation in (3.13) implies the relation Vn = Un that is valid in the GR ACDM context, 
and therefore onr eqnations rednce to those in [22], as they shonld. 

A key pnrpose of the post-Friedmann approach is that it allows ns to derive additional 
eqnations, snch as the eqnation for the vector potential. We now apply the same expansion 
to the Oi Einstein eqnation in order to derive this eqnation. The Oi Einstein eqnation in f{R) 
gravity is given by 


Slhi (1 + faiR) + S/r) -^{r + 6R + f{R) + 6f{R)) 

+goia{fR{R) + 6Jr) - VoViSfR = -^pa(l + S)vi. (3.14) 
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As before, expanding ont the terms and then keeping the leading order only leads ns to 
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As in GR, the scalar part of this final equation is rednndant: If we take the divergence of 
this eqnation then we get an eqnivalent eqnation to the time derivative of eqnation (3.9), 
nsing energy-momentnm conservation to snbstitnte for 6. Note that the continnity and Enler 
eqnations for the matter are the same as for a ACDM cosmology, see [22] for these eqnations 
in the post-Friedmann approach. 

Taking the vector (divergence free) part of the final eqnation in (3.15) leads to the 
following eqnation for the vector potential at leading order. 


= - 


IGvrGpa 
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This is the same as the eqnivalent eqnation in a GR-I-ACDM cosmology, which is easily 
nnderstood as follows. At leading order in this weak-field regime, any additional terms that 
contain the Jr will be linear in the field. Since the new field is a scalar, no combination of 
derivatives of this field can create a divergenceless qnantity withont mnltiplying the field by 
itself or a metric potential, which cannot occnr at order c~^. Thns, the modified Einstein 
eqnations can only modify the scalar part of the Oi eqnation at leading order and not the 
vector part. We expect that this resnlt conld hold for all modified theories of gravity that 
only contain scalar additional degrees of freedom. Nonetheless, the momentnm field itself will 
evolve differently, so the sonrce term for the vector will have a different nnmerical valne. 


4 Vector potential from f{R) N-body simulations 

Using eqnation (3.16), we can extract the vector potential from f{R) N-body simnlations. 
We will follow the extraction process as ontlined in [3, 31] for the momentnm field method. 
This consists of nsing the Delanney Tesselation Field Estimator (DTFE) code [6, 27, 32] to 
extract the sonrce term for the vector potential, the momentnm field. After extracting the 
momentnm field from the simnlations, the scalar part is snbtracted in fonrier space, 

[{i + i)v]|„ = (i + «)v-k2G(l±Ml, (41) 

The power spectrnm of the vector potential is then given by 

ni»(t) = (l^^^) (4.2) 

The momentnm field can also be extracted (and the vector potential calcnlated) nsing a 
Glond-in-Gells code [11] for comparison. As described in [31], there is also an additional way 
to extract the power spectrnm of the vector potential nsing the DTFE code, which revolves 
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around taking the curl of equation (3.16). This method has also been applied to the f{R) 
snapshots with reasonable agreement on small scales, but worse agreement on larger scales. 
The results from these methods are detailed in appendix A. For this work, the simulations 
were run with 1024^ particles in a 250/i“^Mpc box. The initial conditions were created at 
redshift 49 using the Zel’dovich approximation [33]. The particles were then evolved to z = 0 
using the ECOSMOG code [19] (see also [15, 18] for more details), which was based on 
the RAMSES code [29]. The particles were evolved under the influence of f{R) gravity for 
two different values of the parameter, |/ro| = 1.289 x 10“® and |/rq| = 1.289 x 10“®. 
The same initial distribution was also evolved under standard GR for comparison. The 
cosmological parameters used for the simulation are as follows: = 0.267, Oa = 0.734, 

Hq = 71kms“^Mpc“^, ag = 0.8 and Ug = 0.958. From [3, 31], we know that the vector 
potential requires a good mass resolution for the calculation of the vector potential. The 
resolution of the simulations used here is not quite in the ideal parameter range, so it is 
possible that there may be a small spurious contribution to the vector potential power spectra 
extracted here. We aim to ameliorate this effect by examining the ratio of the value in the 
f{R) and GR simulations. We have examined the output of the simulations at both redshift 
zero and redshift one, and checked that the results from our GR simulation agree with the 
results obtained in [31]. 

4.1 Results 

Our results for the vector potential are shown in figures 1 and 2. The first of these shows 
the dimensionless vector potential power spectrum at redshift zero, for GR (black, dashed), 
l/i^ol = 1.289 X 10“® (red) and |/ro| = 1.289 x 10“® (blue). Figure 2 shows the ratio of the 
vector potentials extracted at redshift one and redshift zero for the three models. Although 
the vector potential is growing on non-linear scales for all three models, on linear scales the 
vector potential is actually decreasing, similarly to the linear theory scalar potential. Note 
that in all cases the vector potential is much smaller than the scalar gravitational potential 
and shows a similar scale dependence (see [31]). 

In figures 3 and 4 we show the ratio of the vector potential power spectrum in f{R) 
gravity (for I/rqI = 1.289 x 10“® and |/roI = 1-289 x 10“® respectively) to that in GR. In 
these plots, the black curve shows the ratio at redshift zero and the blue curve shows the 
ratio at redshift one. We can see that the vector potential power spectrum in f{R) gravity 
is sensitive to the value of , with the difference compared to GR changing from of order a 
few percent to of order a few tens of percent between the two different values of //?„. In both 
cases, there is a smaller difference compared to GR at earlier times and on larger scales, in 
line with the behaviour of the density and velocity fields in f{R) gravity [18, 19]. 

Paralleling the GR case examined in [3, 30, 31], there are several immediate consequences 
that can be drawn from the value of the vector potential in f{R) gravity. As in GR, the vector 
potential is the leading order correction to the equations valid in the Newtonian-like regime, 
so the small magnitude of this quantity suggests that the approximations used in deriving the 
equations used in f{R) N-body simulations are indeed satisfied to a high degree of accuracy. 
Furthermore, as noted in [30], the equation for the deflection angle in section IIIA is valid for 
any metric theory of gravity, including f{R) gravity. Thus, similarly to GR [3, 30], the vector 
potential in f{R) gravity is unlikely to have a significant effect on weak-lensing observations. 
This also means that the standard ray-tracing approach to weak leasing should be valid in 
f{R) gravity. Nonetheless, the larger amplitude of this vector potential in f{R) gravity on 
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Figure 1. The dimensionless power spectrum of the vector potential extracted from the GR and f{R) 
N-body simulations at redshift zero. The red (solid) curve is for f{R) gravity for | = 1.289 x 10“^ 
and the blue (dotted) curve is for \fRg \ = 1.289 x 10“®. The dot-dashed black curve shows the power 
spectrum in GR. 



k/h Mpc"‘ 


Figure 2. The ratio of the vector potential power spectra at redshift zero and redshift one. The 
red (solid) curve is for f{R) gravity for \fRg\ = 1.289 x 10“® and the blue (dotted) curve is for 
|/j?ol = 1-289 X 10“®. The dot-dashed black curve shows the ratio for GR. 


smaller scales means that, if a way can be found to observe this quantity, it can be used as 
an additional test of modified gravity theories on non-linear scales. 

Although little work has gone into examining f{R) gravity beyond linear perturbation 
theory, the vector potential examined here could in principle also be calculated using second 
order perturbation theory, as done for GR in [14, 20]. Perturbation theory would be expected 
to break down at most of the scales considered here, and is wrong by nearly two orders 
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Figure 3. The ratio of the vector potential power spectrum in f{R) gravity to that in GR, for 
I/roI = 1-289 X 10“^. The blue (dotted) curve shows the ratio at redshift one, and the black curve 
shows the ratio at redshift zero. 



Figure 4. The ratio of the vector potential power spectrum in f{R) gravity to that in GR, for 
|/i?ol = 1.289 X 10“®. The blue (dotted) curve shows the ratio at redshift one, and the black curve 
shows the ratio at redshift zero. 


of magnitude for GR on the smallest scales considered in [3, 31]. However, on the largest 
scales looked at here, we are close to the linear regime, so the vector potential would be 
expected to agree reasonably well with the prediction from perturbation theory. In addition, 
the vorticity could also be examined analytically, which is one of the contributions to the 
alternative curl method of extracting the vector potential (see appendix A) and which has not 
yet been examined in f{R) gravity. It would be interesting to perform detailed perturbative 
calculations in f{R) gravity in order to examine these phenomena and compare them to the 
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results presented here. We leave this for future work. 


5 Conclusions 

In this paper we have applied the post-Friedmann approach to f{R) gravity, and in doing so 
have provided a coherent derivation of the equations used in f{R) N-body simulations. In 
additon, we have derived the equation for the hrst correction to the leading order equations, 
the vector potential. The equation for this vector potential is the same as in GR, a result that 
may hold for other modihed gravity theories that only contain an additional scalar degree of 
freedom. 

We have used this equation to extract the power spectrum of the vector potential from 
f{R) N-body simulations, following the procedure used for a GR-I-ACDM cosmology in [3, 31]. 
The result of this analysis is shown in hgures 1-2 and hgures 3-4. We found that the vector 
potential can be close to 50% larger in the f{R) simulations for |/ro| = 1.289 x 10“^, although 
this is reduced at earlier times, on larger scales, or for a smaller value of |/rq|. The relatively 
low magnitude of the vector potential suggests that the approximations used in deriving the 
equations for the N-body simulations are accurate and that the vector potential is unlikely 
to be found through weak-lensing surveys. 

The results of this analysis, and that of [3, 31], suggests that even on scales where 
the density contrast is not constrained to be small, gravitational phenomena can be entirely 
described by the two scalar potentials in the metric and their leading order relationship to 
the matter helds. This inference should be checked for further alternative gravity theories 
but, if found to continue to hold, could have important consequences. In particular, although 
much work has been put into model independent parameterisations of modihed gravity, to 
date these parameterisations are perturbative and thus valid only on linear scales. If our 
result holds for other theories of modihed gravity, this could be used as the basis for a model 
independent parameterisation of modihed gravity on non-linear scales. This would be of great 
utility for the analysis of future surveys including the Euclid satellite. 

Whilst this manuscript was being prepared, [7] appeared on the arXiv, which contains a 
similar c~^ expansion applied to f{R) gravity theories. The Hu-Sawicki model studied here 
is in their Glass I theories, for which their results seem equivalent to those obtained for the 
c~^ expansion here. 
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A Additional methods for extracting the vector potential 

As examined in [31], there are two methods for extracting the vector potential from N-body 
simulations. Here we present results from several of those for comparison with our result in 
section 4. Firstly, it is possible to extract the momentum held using a Cloud-in-Cells [11] code 
rather than the DTFE code. Figure 5 shows the ratio between the vector potential calculated 
using these two methods. In this hgure, the black curve shows the ratio for GR and the blue 
curve shows the ratio for f{R), although the two curves are almost indistinguishable. The 
agreement between the two methods is good for most of the range of scales examined here, 
although it becomes worse for the smallest scales. This is expected due to the different window 


- 13 


1.4 



k/h Mpc"' 


Figure 5. The vector potential calculated from the Cloud-in-Cells momentum field divided by that 
calculated from the DTFE momentum field. The black curve is for GR and the blue (dotted) curve 
is for f{R). 

functions that come into the two methods: The DTFE and CiC methods have a different 
method of smoothing the particles in the snapshot onto a regular grid, so the extracted field 
is a convolution of the field with the respective window function. This begins to have a non- 
negligible effect on the extracted field as we move closer to the Nyquist frequency of the grid, 
which is what is being seen here. 

In addition, rather than extracting the full momentum field and then calculating the 
vector part, we can take the curl of equation (3.16) as done for GR in [31]. The source term 
then breaks down into three components, VJ x v, 5V x v and V x v. The advantage of 
splitting up the source like this for GR was two-fold: Firstly, it allowed us to compare the 
vorticity we are extracting with the literature, as there were no similar results to compare 
to when the work of [31] was carried out. Furthermore, it allowed us to show that the 
non-linear quantities were contributing far more significantly than the vorticity, which only 
contributes linearly to the source term. Although these advantages are negated for f{R) 
gravity due to the lack of studies of the vorticity in the literature, we can still extract the 
vector potential using this method. As shown in [31], this method does not agree exactly 
with the momentum field method used in the main body of the paper, with the agreement 
being worse on larger scales. In figure 6, we show the difference between the two methods for 
GR (black) and f{R) (blue); the agreement for GR agrees with that expected from [31], with 
the f{R) agreement being even worse on large scales. We are unsure why the two methods 
do not agree better, however we note that the difference between the methods is insufficient 
to affect our conclusions regarding the validity of the approximation or the observability of 
the vector potential. 

In figure 7, we show the analogue of figure 3 for this “curl” method of extracting the 
vetor potential, so the black curve is the ratio between the f{R) (|/roI = 1.289 x 10“^) and 
GR power spectra at redshift z = 0 and the blue curve is the ratio at redshift z = 1. For most 
of the range the result is similar to that obtained in the main body of the paper however, 
at large scales, the increased difference between the two methods for f{R) gravity creates an 
increased ratio between the f{R) and GR result. It seems unlikely that this result is physical. 
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Figure 6. The ratio of the vector potentials calculated using the momentum-field (MF) method and 
the curl method, for both GR (black) and f{R) (blue). 



Figure 7. The ratio between the f{R) (|/i?o| = 1.289 x 10“®) and GR vector potential power spectra 
using the curl method. The black curve shows the ratio at redshift zero, and the blue (dotted) curve 
shows the ratio at redshift one. 


since f{R) quantities typically return to those of GR on these scales [18, 19]. 








